Scenario of strongly non-equilibrium Bose-Einstein condensation 
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Large scale numerical simulations of the Gross-Pitaevskii equation are used to elucidate the self- 
evolution of a Bose gas from a strongly non-equilibrium initial state. The stages of the process 
confirm and refine the theoretical scenario of Bose-Einstein condensation developed by Svistunov, 
Kagan, and Shlyapnikov ||l|, H, H]: the system evolves from the regime of weak turbulence to superfluid 
turbulence via states of strong turbulence in the long-wavelength region of energy space. 
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I. INTRODUCTION 
A. Statement of the problem 

The experimental realization of Bose-Einstein conden- 
sates (BEC) in dilute alkali and hydrogen gases Q and 
more recently in a gas of metastable helium [^ has stimu- 
lated a great interest in the dynamics of BEC. In the case 
of a pure condensate, both equilibrium and dynamical 
properties of the system can be described by the Gross- 
Pitaevskii equation (GPE)[g| (in nonlinear physics this 
equation is known as defocusing nonlinear Schrodinger 
equation). The GPE has been remarkably successful in 
predicting the condensate shape in an external potential, 
the dynamics of the expanding condensate cloud, the mo- 
tion of quantized vortices; it is also a popular qualitative 
model of superfluid helium. 

An important and often overlooked feature of the GPE 
is that it gives an accurate microscopic description of the 
formation of BEC from the strongly degenerate gas of 
weakly interacting bosons IQ, ||- By large scale numeri- 
cal simulations of the GPE it is possible, in principle, to 
reveal all the stages of this evolution from weak turbu- 
lence to superfluid turbulence with a tangle of quantized 
vortices as was argued by Svistunov, Kagan, and Shlyap- 
nikov P, 0, 0] (for a brief review see Ref. pi) . This task has 
up to now remained unfulfilled, though some important 
steps in this direction have been made in Refs. ||l^, pT| . 
We would also like to mention the description of the 
equilibrium fluctuations of the condensate and highly oc- 
cupied non- condensate modes using the time-dependent 
GPE (H Ol- 

The goal of this paper is to obtain the conclusive de- 
scription of the process of strongly non-equilibrium BEC 
formation in a macroscopically large uniform weakly in- 
teracting Bose gas using the GPE. We are especially in- 
terested in tracing the development of the so-called co- 
herent regime [Q, ^, ||, ^ at a certain stage of evolu- 
tion. According to the theoretical predictions fi 0], this 
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regime sets in after the breakdown of the regime of weak 
turbulence in a low-energy region of wavenumber space. 
It corresponds to the formation of the superfluid short- 
range order which is the state of superfluid turbulence 
with quasi-condensate local correlation properties. 

The notions of weak turbulence and superfluid turbu- 
lence are crucial to our understanding of ordering kinet- 
ics. In the regime of weak turbulence (for an introduction 
to the weak turbulence theory for the GPE see Ref. [u] ), 
the "single-particle" modes of the field are almost in- 
dependent due to weak nonlinearity of the system. The 
smallness of correlations between harmonics in the regime 
of weak turbulence implies the absence of any order. On 
the other hand, the regime of superfluid turbulence (for 
an introduction, see Ref. Iq) is the regime of strong coher- 
ence where the local correlation properties correspond to 
the superfluid state, but the long-range order is absent 
because of the presence of a chaotic vortex tangle and 
non-equilibrium long- wave phonons H . In the case of a 
weakly interacting, gas the local superfluid order is syn- 
onymous to the existence of quasi-condensate correlation 
properties |g]. In a macroscopically large system, the 
cross-over from weak turbulence to superfluid turbulence 
is a key ordering process. Indeed, in the regime of weak 
turbulence there is no order at all, while in the regime 
of superfluid turbulence the (local) superfluid order has 
already been formed. Meanwhile, a rigorous theoretical 
as well as numerical or experimental studies of this stage 
of evolution have been lacking. The general conclusions 
concerning this stage ||, || were made on the basis of 
qualitative analysis that naturally contained ad hoc el- 
ements. The difficulty with an accurate analysis of the 
transition from weak turbulence to superfluid turbulence 
comes from the fact that the evolution between these two 
qualitatively different states takes place in the regime of 
strong turbulence which is hardly amenable to analytical 
treatment. Large computational resources are necessary 
for a numerical analysis of this stage since the problem in- 
volves significantly different length scales and, therefore, 
requires high spatial resolution. 

In the present paper we demonstrate that this prob- 
lem can be unambiguously solved with a powerful enough 
computer. Our numerics clearly reveal the dramatic pro- 
cess of transformation from weak turbulence to superfluid 



turbulence and, thus, fills in a serious gap in rigorous the- 
oretic description of the strongly non-equilibrium BEC 
formation kinetics in a macroscopic system. 

The paper is organized as follows. In Sec. I B we discuss 
the relevance of the time-dependent GPE to the descrip- 
tion of the BEC formation kin etics and its relation to the 
other formalisms. In Sec. I C we render some important 



details of the evolution scenario that we are going to ob- 
serve. In Sec. p we describe our numerical procedure. 
In Sec. [II we present the results of our simulations. In 
Sec. IV we conclude with outlining the observed evolu- 



tion scenario and making a comment on the case of a 
confined gas. 



Time-dependent Gross-Pitaevskii equation and 
BEC formation kinetics 



In this Section we will discuss the question of appli- 
cability of the GPE to the BEC formation kinetics, and 
its connections to the other — full-quantum — treatments. 
These discussion is especially relevant in the wake of a re- 
cent controversy on the applicability of the classical-field 
description to a non-condensed bosonic field. 

A general analysis of the kinetics of a weakly interact- 
ing bosonic field was performed in Ref. Eq. In terms of 
the coherent-state formalism, it was demonstrated that 
if the occupation numbers are large and somewhat un- 
certain (with the absolute value of the uncertainty be- 
ing much larger than unity and with a relative value of 
the uncertainty being arbitrarily small), then the system 
evolves as an ensemble of classical fields with correspond- 
ing classical-field action. (For an elementary demonstra- 
tion of this fact for a weakly interacting Bose gas and 
especially for the discussion of the structure of the ini- 
tial state see Ref. ||.) This has a direct analogy with the 
electro-magnetic field: (i) the density matrix of a com- 
pletely disordered weakly interacting Bose gas with large 
and somewhat uncertain occupation numbers is almost 
diagonal in the coherent-state representation, so that the 
initial state can be viewed as a mixture or statistical en- 
semble of coherent states; (ii) to the leading order each 
coherent state evolves along its classical trajectory which 
in our case is given by the GPE 
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where -0 is the complex- valued classical field that specifies 
the index of the coherent state, m is the mass of the bo- 
son, U = A-kU a/m is the strength of the (5-function inter- 
action (pseudo-)potential, and a is the scattering length. 
[Note that in a strongly interacting system it is impossi- 
ble to divide single-particle modes into highly occupied 
and practically empty ones, so the requirement of weak 
interactions is essential here. In a strongly interacting 
system, there are always quantum modes with occupa- 
tion numbers of order unity that are coupled to the rest 
of the system.] Therefore, the behavior of the quantum 



field is equivalent to that of an ensemble of classical mat- 
ter fields. 

It is important to emphasize that in the context of 
the strongly non-equilibrium BEC formation kinetics the 
condition of large occupation numbers is self-consistent: 
the evolution leads to an explosive increase of occupa- 
tion numbers in the low-energy region of wavenumber 
space Q where the ordering process takes place. Even 
if the occupation numbers are of order unity in the ini- 
tial state, so that the classical matter field description is 
not yet applicable, the evolution that can be described at 
this stage by the standard Boltzmann quantum kinetic 
equation inevitably results in the appearance of large oc- 
cupation numbers in the low-energy region of the particle 
distribution (see, e.g., Ref. [ff^). The blow-up scenario 
[|l[ indicates that only low-energy part of the field is ini- 
tially involved in the process. Therefore, one can switch 
from the kinetic equation to the matter field descrip- 
tion for the long-wavelength component of the field at 
a certain moment in the evolution when the occupation 
numbers become appropriately large. As the time scale 
of the formation of the local quasi-condensate correla- 
tions is much smaller than any other characteristic time 
scale of evolution [0, the cut-off of the high-frequency 
modes, associated with the matter field description, is 
not important. By the time the interactions (particle ex- 
change) between the high- and low-frequency modes be- 
came significant, the local superfluid order had already 
been developed. The order of interaction wavelengths 
is of typical thermal deBroglier wavelength and, there- 
fore, these interactions are essentially local with respect 
to the quasi-condensate and can be described in terms of 
the kinetic equation|]ll O. 

The thesis of the applicability of the matter field de- 
scription at large occupation numbers was justified by the 
analysis of Ref. |l^. Later Stoof questioned the validity of 
this thesis by introducing the concept of "quantum nu- 
cleation" of condensate as a result of an essentially quan- 
tum instability pq ] ; the path- integral version of Keldysh 
formalism was used to substantiate this concept. For 
a criticism of the concept of "quantum nucleation" see 
Refs. I, |l|. 

It is important to emphasize, however, that the path- 
integral approach developed in ||l^ appears to be the 
most fundamental, powerful, and universal way of de- 
riving the basic equations for the dynamics of a weakly 
interacting Bose gas. In particular, we believe that the 
demonstration of the applicability of the time-dependent 
GPE to the description of highly occupied single-particle 
modes of a non-condensed gas within this formalism 
would be the most natural since the effective action for 
the bosonic field is simply the classical-field action of the 
GPE. Basically, one simply has to make sure that for 
the modes with large and somewhat uncertain occupa- 
tion numbers the main contribution to the path integral 
comes from a close vicinity of the classical trajectories 
with the quantum corrections being relevant only on large 
enough times of evolution. 



An interesting all-quantum description of the BEC ki- 
netics was implemented in Ref. |ll|. This technique is 
based on associating the quantum-field density matrix in 
coherent-state representation with a correlator of a pair 
of classical fields which evolution is governed by a sys- 
tem of two coupled nonlinear equations with stochastic 
terms. Using this method the authors performed a nu- 
merical simulation of the BEC formation in a trapped 
gas of a moderate size. We believe (in particular, in 
view of general results of Refs. H, |lj) that this approach 
might be further developed analytically to demonstrate 
explicit overlapping with the other treatments and with 
the time-dependent GPE. Indeed, the form of the system 
of two coupled equations of Ref. O is reminiscent of that 
of the GPE. This suggests that under the condition of 
large occupation numbers the system can be decoupled 
leading to the GPE for the diagonal part of the density 
matrix with relative smallncss of the non-diagonal terms. 
If the standard Boltzmann equation is applicable, so that 
the system can be viewed as an ensemble of weakly cou- 
pled elementary modes [Q, it is natural to expect that 
the equations of Ref. O should lead to the kinetic equa- 
tion. A natural way for deriving kinetic equation from 
the dynamic equations of Ref. U^ is to utilize the standard 
formalism of the weak turbulence theory. In the case 
of the GPE, the weak turbulence approximation leads 
to the quantum-field Boltzmann kinetic equation with- 
out spontaneous scattering processes (see, e.g., Ref. [|l4[ ; 
note also that it is the simplest way to make sure that 
the GPE is immediately applicable once the occupation 
numbers are large). It is natural to expect that in the 
full-quantum treatment of Ref. |ll| the weak-turbulence 
procedure over the dynamic equations would result in 
the complete quantum-field Boltzmann kinetic equation 
with the spontaneous processes retained. Unfortunately, 
we are not aware of such investigations of the equations 
of Ref. ^ that might be very instructive for the general 
understanding of the dynamics of a weakly interacting 
Bose gas. 



C. Initial state and evolution scenario 

In what follows we consider the evolution of Eq. (||) 
starting with a strongly non-equilibrium initial condition 



V'(r, t = 0) = 2_] flk exp(ikr). 



(2) 



where the phases of the complex amplitudes Ok are dis- 
tributed randomly. Such an initial condition follows from 
the microscopic quantum- mechanical analysis of the state 
of a weakly interacting Bose gas in the kinetic regime 
Q. Theoretical investigations of the relaxation of such 
an initial state towards the equilibrium configuration 
were performed by Svistunov, Kagan, and Shlyapnikov 
[Q, ||, ^. The analysis revealed a number of stages in 
the evolution. Initially the system is in the weak tur- 
bulence regime and thus can be described by Boltzmann 



kinetic equation. The kinetic equation is obtained as 
the random-phase approximation of Eq. (|l|) for occupa- 
tion numbers rtk defined by (akO^,) ~ Jik'^k'^k'- Alterna- 
tively, the weak-turbulence kinetic equation follows from 
the general quantum Boltzmann kinetic equation if one 
neglects spontaneous scattering as compared with stim- 
ulated scattering (because of large occupation numbers) . 
Svistunov [1| and later Semikoz and Tkachev [|l^ con- 
sidered the self-similar solution of the Boltzmann kinetic 
equation: 



n,{t) = Aeo"(t)/(e/eo), t < t, 

eo(t) = i?(t, -t)i/2("-i), 

f{x) -> a;~" at x ^ cx), /(O) : 



(3) 

(4) 
(5) 



where e = h k"^ /2'm. The dimensional constants A and B 
relate to each other by (a - l)m^U^A'^ = Xir^h'^ B^^°'-^\ 
where the parameters a and A were determined by nu- 
merical analysis as a « 1.24 Q and A « 1.2 ||l|. The 
form of the function / was also determined numerically 
in Ref. [Q. The solution (|^)-(|^) has only one free pa- 
rameter, say A, that depends on conditions in the "pre- 
historic" evolution. By "pre-historic" evolution we mean 
any sort of non- universal dynamics preceding the appear- 
ance of self-similarity. This dynamics is sensitive to the 
details of the initial condition or/and cooling mechanism 
as well as to the spontaneous-scattering terms in the ki- 
netic equation which cannot be neglected until the occu- 
pation numbers are large enough. When the self-similar 
regime sets in at a certain step of evolution all the par- 
ticular details of the previous evolution are absorbed in 
the single parameter A. 

The self-similar solution (§)-(^ describes a wave in 
energy space propagating from high to lower energies. 
The energy eo{t) defines the "head" of the wave. The 
wave propagates in a blow-up fashion: £o{t) —>■ and 
n^g (t) — > cxD as t — > i* . In reality the validity of the 
kinetic equations associated with the random phase ap- 
proximation breaks down shortly before the blow-up time 
t,. This moment marks the beginning of a qualita- 
tively different stage in the evolution the coherent regime: 
strong turbulence evolves into a quasi-condensate state. 
In the coherent regime the phases of the complex ampli- 
tudes Ok of the field ip become strongly correlated and the 
periods of their oscillations arc then comparable with the 
evolution times of the occupation numbers. The forma- 
tion of the quasi-condensate is manifested by the appear- 
ance of a well-defined tangle of quantized vortices and, 
therefore, by the beginning of the final stage of the evo- 
lution: superfluid turbulence. In this regime the vortex 
tangle starts to relax over macroscopically large times. 



II. NUMERIC PROCEDURE 
A. Finite-difference scheme 

We performed a large scale numerical integration of a 
dimensionlcss form of the GPE 



2i^ = VV+IV'P^, 



(6) 



starting with a strongly non-equilibrium initial condition. 
Our calculations were done in a periodic box N^, with 
iV = 256, using a fourth-order (with respect to the spa- 
tial variables) finite-difference scheme. The scheme cor- 
responds to the Hamiltonian system in the discrete vari- 
ables ■tpijk'- 



.dip, 



'ijk 



at 



dH 



(7) 



where 



H = 



Ij: 



t 



ijk 



[—j2{''Pi+2,j,k — fpi-2,j,k 



ijk 



+ V-'i,j+2,fc —1piJ~2,k + 'ipi,j,k+2 — 1pi,j,k-2) 
+ -g{ '4'i+l,j,k —1pi-l,j,k + iJi,j+l,k — 1pi,j-l,k 
+ 'ipi,j,k+l -Ipi.j.k-l)]- ^I'ipijkl'^ (8) 

[in the numerics we set the space step in each direction 
of the grid as dx — dy ^ dz = 1]. 

Equation (R) conserves the energy H and the total par- 
ticle number 2^4,;;, IV'ijfeP exactly. In time stepping, the 
leap-frog scheme was implemented 






t 



ijk 



2dt 



dH 



dt 



ijk 



(9) 



with dt = 0.03. To prevent the even-odd instability of 
the leap-frog iterations, we introduce the backward Euler 
step 



>Sr 



^. 



ijk 



dH \"+i 



dt 






(10) 



every 10"* time steps. The leap-frog scheme is nondissi- 
pative, so the only loss of energy and of the total particle 
number occurs during the backward Euler step and, since 
we take this step very rarely, these losses are insignificant. 
The code was tested against known solutions of the 
GPE: vortex rings and rarefaction pulses [ pi| . The sim- 
ulations were performed on a Sun Enterprise 450 Server 
and took about three months to complete for the main 
set of calculations discussed below. 



B. Initial condition 

To eliminate the computationally expensive (and the 
least physically interesting) transient regime, we started 



directly from the self-similar solution Eq. (||) with (||)-(||), 
so that: 



flk = \/Ck"o/(e/eo) exp[i0k] 



(11) 



where ^k and (/)k are random numbers. [Note that in the 
simulations the momentum k is the momentum of the 
lattice Fourier transform.] The phase 0k is uniformly dis- 
tributed on [0, 27r] in accordance with the basic statement 
of the theory of weak turbulence and with the explicit mi- 
croscopic analysis of corresponding quantum field states 
[p| . The choice of ^k is rather arbitrary, we only fix its 
mean value to be equal to unity, introducing therefore 
the parameter rig. The weak-turbulence evolution is in- 
variant to the details of the statistics of |ak|. By the 
time the system enters the regime of strong turbulence, 
the proper statistics is established automatically since 
each harmonic participates in a large number of scatter- 
ing events. We tried different distributions for ^k and saw 
no systematic difference in the evolution picture. The 
main set of our simulations was done with the distribu- 
tion function t«(^k) = exp (— ^k) (heuristically suggested 
by equilibrium Hibbs statistics of harmonics in the non- 
interacting model). 

When choosing the parameters of the initial condition 
(0) specified by the complex Fourier amplitudes (pT|), we 
have to take eg small enough to be free from systematic 
error of large finite-differences. On the other hand, taking 
Eq too small reduces the physical size of the system. Let 
us define one period of the amplitude oscillation as tp — 
27r/eo and the number of periods before the blow-up as 
P = t^,/tp. When choosing the value of no in combination 
with eo, we would like to avoid having P too small when 
the time scale of the kinetic regime becomes too short, 
or having P too large, when the finite-size effects (the 
discreteness of the k) dominate the calculation. Given 
the maximal available grid size N — 256, we found that 
it is optimal to take no = 15 and eo — 1/18, so that 
tp « 113, t^ = 4A7rV(a - l)egno « 893, and P « 8. 



III. DATA PROCESSING AND RESULTS 

The instantaneous values of the occupation numbers 
"■k(i) = |ak(OI^ ^''"S extremely 'noisy' functions of time. 
To be able to draw some quantitative comparisons and 
conclusions we need either to perform some averaging or 
to deal with some coarse-grained self-averaging charac- 
teristics of the particle distribution. Taking the second 
option, we introduce shells in momentum space. By the 
i-th shell (i = 1,2,3,.. .) we understand the set of mo- 
menta satisfying the condition i — 1 < log2(fc/27r) < i. 
The idea behind this definition is that each shell rep- 
resents some typical momentum (wavelength) scale and 
thus allows to introduce a coarse-grained characteristic of 
the occupation numbers corresponding to a given scale. 
Namely, for each shell i we introduce the mean occupa- 
tion number rjiit) =^ n\^li:,)/Mi, where Mi is the 
number of harmonics in the i-th shell. The harmonic 
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FIG. 1: The time evolution of r]i{t — lOOj) in the weak tur- 
bulence regime for j — 0, ... ,6. In the insets we show the 
theoretical self-similar solution (|3)-(p|) (solid line) and the 
solution obtained through the numerical integration of (n) 
(dashed line) for the shells i = 1 (a) and i — 2 (b). 



k = 0, which plays a special role (at the very end of the 
evolution), is not assigned to any shell. 

Another instructive coarse-grained characteristic of the 
particle distribution is the integral distribution function 
Fk — X]fc'<fc"'k' which shows how many particles have 
momenta not exceeding k. We use function Fk to keep 
track of the formation of the quasi-condensate and to de- 
termine wavenumber span of the above-the-condensate 
particles. This information is used, in particular, for fil- 
tering out the high-frequency harmonics in order to inter- 
pret the results of our numerical calculations in superfluid 
turbulence regime. 

With the above-introduced quantities we now turn to 
the analysis of the results of our numerical simulations. 
The self-similar character of the evolution is clearly ob- 
served in Fig. 1. Inserts in Fig. 1 give the comparison 
of the theoretical prediction of the evolution of the occu- 
pation number function ne{t) defined by Eq. (g) and the 
evolution of the first and the second shells. The agree- 
ment with the theoretical predictions W is quite good for 
t < 600. After that the numerical solution deviates from 
the self-similar theoretical solution which is the mani- 
festation of the onset of the strong turbulence stage of 
evolution. 

As follows from the dimensional analysis (see, e.g., 
Ref. [pi), the characteristic time, io, and the characteris- 
tic wave vector, kg, at the beginning of the strong turbu- 
lence regime are given by the relations 






to - Cq [h^'^+'/m^U^A^] i/(2"-i) , (12) 



fco - Ci[At/(TO/?i)"+l]l/(2a-l) 



(13) 



where Cq are Ci some dimensionless constants. Our 
numerical results (Fig. 1) indicate that t^, — to ~ 300 
which implies Cq ^ 40. After the formation of the quasi- 
condensate {t > 1000), the distribution of particles ac- 
quires a bimodal shape which is seen in Fig. 2. The 



FIG. 2: Evolution of the integral distribution of particles 
Fk = X/fc'<fc ^k'- Notice the appearance of a 'shoulder' of Fk 
indicative of the quasi-condensate formation. The evolution 
of nk=o is presented in the inset. Note the strong fluctuations 
typical for the evolution of a single harmonic. The fluctua- 
tions are also seen in the graph of the first shell (see insert 
(a) of Fig. 1). 



salient characteristic of the distribution is the shoulder 
which becomes sharper and sharper as the evolution con- 
tinues. Note that by definition of the function Fk, the 
height of the shoulder is equal to the number of the quasi- 
condensate particles. From Fig. 2 we estimate ko ^ 15 
which means that 



Ci - Co - 40 



(14) 



Within the coherent regime the momentum distribu- 
tion of the harmonics yields rather incomplete picture 
of the evolution and it becomes reasonable to follow the 
ordering process in coordinate space. The tracing the 
topological defects in the phase of the long-wavelength 
part of the complex matter field ip is very important 
since the transformation of these defects into a tangle of 
well-separated vortex lines is the most essential feature 
of the superfiuid short-range ordering g. To this end 
we first filter out the high-frequency harmonics by per- 
forming the transformation Ok — > aknraxjl — fc^/fc^,0}, 
where kc is a cut-off wave number. When the function Fk 
has a pronounced quasi-condensate shoulder, the natural 
choice is to take kc somewhat larger than the momen- 
tum of the shoulder in order to remove the above-the- 
condensate part of the field ijj. In the regime of weak 
turbulence, when there is no quasi-condensate, the pro- 
cedure of filtering is ambiguous: the distribution is not 
bimodal, so there is no special low momentum k^, also, 
the structure of the defects in the filtered field essentially 
depends on the cut-off parameter and, thus, has no phys- 
ical meaning. 

The results of visualizing the topological defects are 
presented in Fig. 3. The formation of a tangle of well- 
separated vortices and the decay of superfluid turbulence 
are clearly seen. This is the key point of our simulation. 
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FIG. 3: Evolution of topological defects in the phase of the long-wavelength part tp of the field ^ in the computational box 
256^. The defects are visualized by isosurfaces lipf = 0.05(|?/'p). High-frequency spatial waves are suppressed by the factor 
max{l — fe^/fcc, 0}, where the cut-off wave number is chosen according to the phenomenological formula kc — 9 — t/1000. 



To the best of our knowledge, this is the first unambigu- 
ous demonstration of the formation of the state of super- 
fluid turbulence in the course of self-evolution of weakly 
interacting Bose gas. This result forms a solid basis for 
the analysis of the further stages of long-range ordering in 
terms of well-developed theory of supcrfluid turbulence 
that was performed in Rcf. g (see also Rcf. ^). 

The characteristic time of the evolution of the vortex 
tangle depends on the typical interline spacing, R, as 
R'^/\n{R/ao), where ao is the vortex core size (see, e.g., 
Q). During the final stage of evolution, when R is of 
the order of linear size of the computational box, the 
slowing down of the relaxation process makes numerical 
simulation of the final stage of the vortex tangle decay to 
be enormously expensive in a large computational box. 
For example, according to the above-mentioned estimate 
of the relaxation time, to achieve the complete disap- 
pearance of the vortex tangle in our N = 256 system 
we would need several years. To observe this final stage 
of the vortex tangle decay, we repeated the calculations 
for a smaller computational box with A^ = 128 (reduc- 
ing in this way the computational time by a factor of 
~ 32 — 2^ X 2^); see Fig. 4. Parameters of the initial 
condition are e = 1/2 and no = 27r, so that the number 
of periods before the blow-up is P « 5. A single vortex 
ring remains at i = 4000 as the result of the turbulence 
decay; see Fig. 5(a). 



The above-mentioned filtering method allows us to vi- 
sualize the position of the core of a quantized vortex line, 
but not the actual size of the core since we force the so- 
lution to be represented by a relatively small number 
of harmonics. To get a better representation of an ac- 
tual size of the core as well as to resolve another objects 
of interest — rarefaction pulses |2l|, which are likely to 
appear in the course of transformation of strong turbu- 
lence into superfluid turbulence, we implement a different 
type of filtering based on time averaging. We introduce 
a Gaussian- weighted time average of the field ip: 



^zjkit) 



V',jfc(T)exp[-(r-t)VlOO]dr 



(15) 



The width of the Gaussian kernel in (|T^ is chosen in 
such a way that the (disordered) high-frequency part of 
the field ip is averaged out revealing the strongly corre- 
lated low-frequency part ip- Fig. 5 compares the density 
isosurfaces obtained by two different methods: by high- 
frequency suppression (Fig. 5a) and by time averaging 
(Fig. 5b). In the latter case we reveal the actual shape 
of the vortex core and resolve the rarefaction pulses. 



i = 



t = (m 



i=U0O0 




t=iJltO0 



t = aoao 



t=2iW 




FIG. 4: Evolution of topological defects in the phase of the long-wavelength part tp of the field i/' in the computational box 
128^. The defects are visualized by isosurfaces lipf = 0.05(|?/'p). High-frequency spatial waves are suppressed by the factor 
max{l — k^/k^,0}, where the cut-off wave number is chosen according to the phenomenological formula kc — 9 — i/1000. 



IV. CONCLUSION 



We have performed large scale numerical simulations 
of the process of strongly non-equilibrium Bose- Einstein 
condensation in a uniform weakly interacting Bose gas. 
In the limit of weak interaction under the condition of 
strong enough deviation from equilibrium the key stage of 
ordering dynamics — supcrfluid turbulence formation — is 
universal and corresponds to the process of self-ordering 
of a classical matter field which dynamics is governed 
by the time-dependent Gross-Pitaevskii equation (defo- 
cusing nonlinear Schrodinger equation) . The universality 
implies independence of evolution from the details of ini- 
tial processes such as, for example, cooling mechanism 



and rate as well as from quantum effects such as sponta- 
neous scattering. All the information about the evolution 
preceding the universal stage is absorbed in the single pa- 
rameter A that defines scaling of characteristic time and 
wavenumber, in accordance with Eqs. (|l2|)-(|l4|). 

The most important features of the BEG formation 
scenario observed in our simulation are as follows. The 
low-energy part of the quantum field characterized by 
large occupation numbers and, thus, described by a clas- 
sical complex matter field ip obeying Eq. {^ initially 
evolves in a weak-turbulent self-similar fashion accord- 
ing to Eqs. (^)-(l^). The occupation numbers at small 
energies become progressively larger. At the characteris- 
tic time moment tg, given by Eq. (I2h, close to the formal 
blow-up time i* of the solution (P)-(l^) , the self-similarity 
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FIG. 5: Comparison of two isosurfaces obtained by different 
filtering techniques. The solution at f = 4000 is obtained by 
numerical integration of (n) in the periodic box with A'' = 
128. The isosurface \-ii)\^ = 0.05(ii/;|^> is plotted in Fig. 5(a) 
using the high-frequency filtering with fee = 5. The isosurface 
[tpl = 0.2(|i/)| ) is plotted in Fig. 5(b), where ip is defined by 
Eq.(|l|). 



of the energy distribution breaks down. The distribu- 
tion gradually becomes bimodal; the low-energy quasi- 
condensate part of the field sets to the state of superfluid 
turbulence characterized by a tangle of the vortex lines. 
The further evolution of the quasi-condensate is indepen- 
dent of the rest of the system (apart from a permanent 
flux of the particles into the quasi-condensate) and basi- 
cally is the process of relaxation of superfluid turbulence. 
All vortex lines relax in a macroscopically large time. 

In the present paper, we dealt with the case of macro- 
scopically large uniform system. As far as the case of a 



trapped gas is concerned, the situation becomes sensitive 
to the competition between flnite size and nonlinear ef- 
fects. If nonlinear effects dominate, the basic physics of 
the ordering process is predicted to be analogous to that 
revealed by our simulation p^ . If finite-size effects dom- 
inate (which means that the initial size of the condensate 
is smaller than corresponding healing length, so that, for 
example, vortices cannot arise in principle |23|| ), the or- 
dering kinetics is substantially simplified being reduced 
to the growth of genuine condensate ||2^, |4|. Clearly, 
our numerical approach can be extended to the case of 
a trapped Bose gas by simply including a term with an 
external potential in Eq. (11]). Such a simulation could 
provide a deeper interpretation of the first experiments 
on the kinetics of BEC formation |g^, g6| answering, in 
particular, the question of whether the process involves 
the formation of vortex tangle; and if not, under which 
conditions one may expect formation of superfluid turbu- 
lence (quasi condensate) in a realistic experimental situ- 
ation. 
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